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A simple formula is derived for the maximum squeezing rate which occurs at the initial stages 
of the squeezing process: the rate only depends on the second partial derivatives of a classical 
Hamiltonian. Rules for optimum rotation of the phase space are found to keep the state optimally 
located and oriented for fastest squeezing. These operations transform the phase-space point of 
interest into a saddle point with opposite principal curvatures. Similar results are found for the 
Bloch-sphere phase space and spin squeezing. Application of the general formulas is illustrated by 
several model examples including parametric downconversion, Kerr nonlinearity, Jaynes-Cummings 
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I. INTRODUCTION 


Squeezing is an irreducible resource for quantum in¬ 
formation processing [ij. Suppressing the noise of some 
physical variables in squeezed states has important ap¬ 
plications in quantum metrology [H- Sometimes it is 
stressed that squeezed states are purely quantum me¬ 
chanical states as their Glauber P-representation is non¬ 
positive definite Therefore, to find squeezing proper¬ 
ties of various physical systems, it is very natural to use 
mathematical apparatus of quantum physics. 

Relations between classical and quantum predictions 
for various squeezing processes have been studied, dis¬ 
cussing the similarities and differences (see, e.g., [M)- 
However, there has been no general approach showing 
which particular features of classical systems are respon¬ 
sible for noise suppression that would be analogous to 
squeezing production. Here, a simple formula is derived 
for a rate at which squeezing is generated at initial stages 
of the process: the maximum rate only depends on the 
second partial derivatives of the classical Hamiltonian. 
Also, simple formulas are found to determine rotations of 
the phase space by which one keeps the state optimally 
located and oriented to achieve fastest squeezing rate: 
they only contain first and second partial derivatives of 
the classical Hamiltonian. Although for precise results a 
full quantum calculation is necessary, the classical formu¬ 
las work surprisingly well as long as the uncertainty area 
is not deformed significantly beyond an elliptical shape 
described by a variance matrix. The formulas can be used 
as a simple rule of thumb for squeezing prospectors who 
need a quick orientation in the terrain to decide where to 
start mining their precious resource. 

A simple intuitive picture of “classical squeezing” is as 
follows. Imagine a group of tourists starting their hike in 
a hilly countryside. Each member of the group goes along 
a contour line of constant elevation, having the hill on the 
left and the valley on the right, with a speed proportional 
to the magnitude of the slope. Even though initially the 
group might have a circular form, moving on the uneven 
landscape changes the formation to be stretched in one 


direction and squeezed in another. Here, of course, the 
countryside is a phase space, elevation is the value of 
the Hamiltonian, and the hiking rules are the classical 
Hamilton equations. The group of tourists represents an 
ensemble of classical states, and our task is to infer from 
the local shape of the landscape the rate at which the 
group gets squeezed. 

The paper is organized as follows. In Sec. |lT]the for¬ 
mula for the squeezing rate is derived. In Sec. IHII we 
study the question how the uncertainty area changes ori¬ 
entation, in Sec. m we show how to compensate the 
motion of the uncertainty ellipse to keep the optimum 
squeezing rate, and in Sec. |V] we illustrate the general 
results on several examples. In Sec. EH we show how the 
results can be generalized for spin squeezing and motion 
on the Bloch sphere, and we conclude in Sec. EUl Several 
derivations of technical nature are given in Appendixes. 


II. SQUEEZING RATE 


Consider a classical system described by a Hamiltonian 
H(x,p) where the quantities x and p are rescalled such 
that they have the same dimension. We consider a prob¬ 
ability distribution p{x,p) characterized by a variation 
matrix 


/ (Ax^) {AxAp) 
{AxAp) (Ap^) 


^xx ^xp \ 

lAp ^pp J 


( 1 ) 


for a state centered in (xo,po). The question is what fea¬ 
tures of the Hamiltonian determine the squeezing gener¬ 
ation in the system. 

Assume that at time t = 0 the system is in state (xq + 
Ax,po + Ap). At short time dt the system will be in a 
new state (xq + Ax,po + Ap), where up to the first order 
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Denoting the partial derivatives at (a;o,po) as indexes, 
dH{xo,po)/dx = Hx, etc., we can write for the new cen¬ 
tral positions in the phase space 


FIG. 1: (a) Uncertainty ellipse corresponding to variation 
matrix V with the main axis oriented at 6 is transformed into 
a new ellipse corresponding to V with the main axis oriented 
at 6. (b) The squeezing process described by the matrix S 
of Eq. (I19II corresponds to a rotation of the phase space by 
d, stretching along Ax and squeezing along Ap, and rotation 
of the phase space back by (p — e. Quantities x and p are 
dimensionless here and in the next figures. 


(AxAp) = 0, where we get 


(Ai^) s 

s cr^ (1 + 2Hxp) dt, 

(13) 

(Ap2) s 

s cr^ (1 - 2Hxp) dt, 

(14) 

(AiAp) s 

S <7 (^Hpp Hxx) ^t. 

(15) 


io ~ a;o + Hpdt, 
pQ ~ Po - Hxdt, 


(4) To find the rate of squeezing generation, we express the 

(5) eigenvalues of the new variance matrix V as 


and for the new deviations 


Ax R 

:! Ax + {HxpAx HppAp) dt. 

(6) 

Ap R 

i Ap - {HxpAp + HxxAx) dt. 

(7) 

Assuming (Ax) 

= (Ap) = 0 and expressing the new 

vari- 

ances up to the first order in dt we get 


(Ai^) ^ 

. (Ax^) 



-1-2 (Hxp{Ax^) -1- Hpp(AxAp)) dt, 

(8) 

(Ap2) ^ 

. (Ap2) 



-2 (HxpiAp^) + Hxx(AxAp)) dt. 

(9) 

(AiAp) R: 

i (AxAp) 



+ (HppiAp^) - Hxx{Ax^)) dt. 

(10) 


These results can be described as transformation of the 
variation matrix according to 

V = SVS^, (11) 

where 



_ (Ai^) -b (Ap2) 

2 



((Ax^) — (Ap2))^ -|-4(AxAp)2, 

(16) 

finding 


V± = a^(l± Qdt) 

(17) 

where 


Q = ^{Hpp-Hxx)^ + mlp 

(18) 


is the squeezing rate. Note that this formula is invariant 
with respect to rotations of the phase space. In points 
of zero gradient, = Hp = 0, it has the following 
geometric interpretation. If H is taken in the same units 
as X and p, then the principal curvatures of its graph 

are \[Hxx + Hpp) ± \\J[Hpp - Hxx)^ + ^H^p- Thus, in 
this case, Q is proportional to the difference of principal 
curvatures of the Hamiltonian graph. 


/ 1 + Hxpdt Hppdt \ 

y Hxxdt 1 Hxpdt J 


(12) III. ORIENTATION AND ROTATION OF THE 
SQUEEZING ELLIPSE 


with all terms taken up to the first order in dt. 

The results are simple for initially isotropic and un- As can be seen, in the special case of Hxx = Hpp = 0 
correlated fluctuations, i.e., (Ax^) = (Ap^) = cr^, and the transformation matrix S of Eq. (TT^ is diagonal. 
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If also V is diagonal (i.e., Ax and Ap are uncorre¬ 
lated), the transformation squeezes one of the variables 
and stretches the other with rate Q of Eq. (HI), i.e., 
Q = 2\Hxp\. In a general case, however, the squeezing 
ellipse changes orientation of the main axis, as shown in 
Fig. H. This process can be described as follows (see 
Fig. H): the variation matrix is rotated by <() to a new 
coordinate system where S is diagonal, then squeezing 
and stretching occurs along the new coordinates, and the 
variation matrix is rotated back by a modified angle (j)—e. 
Thus, the transformation matrix S can be written as 

c _ f cos {(j)- e) - sin ((/) - e) \ / 1 -I- ^ 0 \ 

“ sin ((() - e) cos (0 - e) J \ 0 1 - ^ / 

(19) 

sm(p COS0 J ' ' 

Expanding this expression up to the first order in e and 
dt, one finds 


( 1-I-cos 2(() e-I-sin2(/) \ 

^ = • ( 20 ) 
\-e+ ^sin 24 > I - ^ cos2(j) J 


Comparing this with Eq. (HH) one finds 


This yields 


Qcos2(j) = 

2Hxp, 

(21) 

Q sin 2(j) = 

Hpp Hxx 5 

(22) 

2e = 

{Hxx Hpp) dt. 

(23) 

tan2(/) = 

Hpp Hxx 

2Hxp ’ 

(24) 


and assuming that e evolves with time as e = uiydt, one 
gets 


(jJy 




H, 


pp 


(25) 


with Q given by Eq. (fT51) . Eq. tells us what is the 
best orientation of the main axis of the uncertainty ellipse 
to generate squeezing the fastest way, namely, 9 = (J). Eq. 
(1^51) tells us with what rate should one rotate the system 
to keep the uncertainty ellipse optimally oriented, namely 

U! = —LOy. 


IV. COMPENSATION OF MOTION OF THE 
UNCERTAINTY ELLIPSE 


During the evolution the uncertainty ellipse not only 
deforms, but also drifts through the phase space. It is 
convenient to express the motion of the center of the el¬ 
lipse as rotation around a phase space point {xr^pb) with 
angular velocity Wc (see Fig. [21). We show in Appendix 


P 



(xr’Pr) 


X 


FIG. 2: Motion of the uncertainty ellipse as combination of 
rotation of its center around {xr,pr) with rate ujc, and change 
of its orientation with rate uiv 


E that for a phase space point {x,p) the rotation center 
is at {xr,ph) = {x,p) + {Rx,Rp) with 


Hx 

ry _ 

— 

u;c 

r, _ ^P 

lip - 


(26) 

(27) 


with the angular frequency of the motion of the center 
being 


HiH, 


UJc = 


PP 


HpHxx - 2HxHpHxp 


HI + Hj 


(28) 


While the center rotates around {xr,pr) with uJc, an 
uncertainty ellipse with optimally oriented main axis de¬ 
forms and changes orientation with angular velocity cOy 
(see Fig. [2]). These formulas could be useful if one is 
able to construct quadratic Hamiltonians of the form 
H = ^uj[{x — Xc)‘^ + {p—pc)'^] with variable parameters w, 
Xc and Pc- This can be achieved, e.g., in quantum opti¬ 
cal experiments where rotations around the phase-space 
origin corresponds to the accumulation of interferometric 
phase, and rotations around other points can be realized 
by combinations of interferometric phase shifts and dis¬ 
placements realized by mixing the quantum field with a 
strong coherent signal on an unbalanced beam splitter 

El- 

Suppose we initiate the system in a state centered at 
(xqjPo) for which Q reaches the desired value. We want 
to keep the state centered here and also keep the un¬ 
certainty ellipse optimally oriented during the squeezing 
process. To compensate for the motion of the uncertainty 
ellipse center we first add the Hamiltonian i7adi in the 
form 

Hadl=-^U}c[{x-XRf + {p-pRf] (29) 


with (xr,pr) calculated according to the above formu¬ 
las and Eqs. (EH), EH) in (x,p) = (a;o,po). Hamiltonian 
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+ _ffadi has zero gradient so that the uncertainty ellipse 
stays centered at (xo,Po)- If hs main axis is at the begin¬ 
ning optimally oriented, it starts rotation with angular 
velocity uj^ — ujc- To keep the optimal orientation, one 
adds another Hamiltonian i/ad 2 in the form 

- Wc) [(x - Xof + [p-pof] (30) 

which rotates the phase space around (a;o,po) with the 
appropriate frequency. As the result, (xo,po) becomes 
a saddle point with principal curvatures of equal mag¬ 
nitude and opposite signs. The additional Hamiltoni¬ 
ans combine to a single quadratic Hamiltonian iJad = 
ffadi + Had 2 SO that the system evolves under the Hamil¬ 
tonian H + Hg_d with 

Had = - 2 ^v [(a; - Xr)^ + {p- Prf] + const., (31) 
where the center is localized at 
{Xr,Pr) = {xr,Pr) 

+ (l -— ] {xo - xr,po - pr). (32) 

\ / 

As can be checked, the squeezing rate Q is unchanged. 

V. EXAMPLES 
A. Harmonic oscillator 

The Hamiltonian is 

H =hx(j>^ + x^), (33) 

and Eq. (IT^ yields (5 = 0, i.e., the harmonic oscillator 
does not produce squeezing. Angular frequencies of Eq. 
(1^ and (1^ are = Wc = w, i.e., equal to the oscillator 
frequency. The motion of the uncertainty ellipse is shown 
in Fig. m- 

B. Quadratic Hamiltonians, parametric down 
conversion 

Apart from the harmonic oscillator, other Hamiltoni¬ 
ans quadratic in x and p generate squeezing. In quantum 
optics they all have a simple interpretation of a paramet¬ 
ric amplifier [3 [11. 

1. Free particle 

With dimensionless x and p, the Hamiltonian of a free 
particle is 

« = Ap". (M) 



FIG. 3: Transformation of the uncertainty ellipse for the har¬ 
monic oscillator (a), free particle (b), inverted oscillator (c), 
and the xp Hamiltonian (d). 


This means that H^x = Hxp = 0 and Hpp = 1/m which 
leads to the squeezing rate Q = 1/m independent of the 
localization in the phase space. The optimum orienta¬ 
tion of the uncertainty ellipse is 0 = 7r/4 which is being 
rotated with the rate uiy = l/(2m), and lUc = 0, see Fig. 
Wp- This means that to keep the optimum orientation, 
one has to rotate the system phase space with the rate 
-l/(2m). 

The quantum optical interpretation of this Hamilto¬ 
nian is found on assuming x and p to be operators con¬ 
structed as combinations of creation and annihilation op¬ 
erators, namely (assuming m = 1) 

X = —j= (fll -I- a) , (35) 

v2 

P = 

The Hamiltonian is then 



The first term on the right hand side corresponds to a 
parametric down conversion with photons being created 
and destroyed in pairs, whereas the second term corre¬ 
sponds to a harmonic oscillator with frequency ujy = 1/2. 
Using the additional Hamiltonian of Eq. (ED) means just 
removing this second term. Note that the squeezing rate 
Q exactly corresponds to the quantum mechanical result 
discussed, e.g., in Q. 
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8. Inverted oscillator 
The Hamiltonian is 


H = lap" - xa ( 38 ) 

which leads to the squeezing rate Q = 2^, independent 
of the initial state. The optimum orientation is 0 = 7r/4, 
and the ellipse does not rotate, a;„ = 0. The evolution of 
the uncertainty ellipse is shown in Fig. |3j:. 

In terms of quantum optical operators the Hamiltonian 
can be written as 

H = laf - x^ =-IC , (39) 

corresponding to the parametric down conversion. 




PP' 

5 

w 

P 


0 


-5 

• 


■ \n\.\ , , . /, 


FIG. 4: (Color online) Contour lines of the Kerr Hamiltonian 
and evolution of the uncertainty lines, (a) Hamiltonian (1441) 
with X = 1, (b) the same Hamiltonian with the added term 
(153. The initial state is centered at {xo,po) ~ (7,0) and the 
uncertainty lines are plotted in times t = 0, 2 x 10“®, 4 x 10“®, 
and 6 x 10~®. 


3. xp-Hamiltonian 

The Hamiltonian is in the form 

H = Cxp (40) 

which is a classical counterpart of the quantum operator 

H = lc{xp + px) 

= - a^) (41) 

corresponding to the parametric down conversion dis¬ 
cussed in detail in Q. Compared to the preceding two 
cases it has just different phase ratio of the quadratures 
X and p. For the Hamiltonian (I40p the squeezing rate 
is Q = 2<C, the optimum orientation is 0 = 0, and no 
rotation is generated, = 0. The evolution of the un¬ 
certainty ellipse is shown in Fig. [3Ji. 

C. Pendulum 

The Hamiltonian is 

H = -p^ — cosx (42) 

leading to the squeezing rate 

Q = 1 — cos X = 2 sin^ ^, (43) 

which changes continuously between 0 for a; = 0 (i.e., 
like harmonic oscillator near the stable equilibrium) and 
2 for X = TT (i.e., like inverted oscillator near the unstable 
equilibrium). The optimum orientation is 0 = 7r/4 and 
the rotation frequency is = cos^ § changing contin¬ 
uously from 1 near the stable equilibrium to 0 near the 
unstable equilibrium. 


D. Kerr uouliuearity 

Assume the Hamiltonian 

H = xip^+xa^ (44) 

whose quantum counterpart was shown to generate 
squeezing 01I1- Eq. (HU) leads to the squeezing rate 

Q = 8x(p^ + xa, (45) 

i.e., the squeezing rate increases with the oscillator en¬ 
ergy. Eq. (HSl) corresponds to the analytical result for 
a quantum Kerr oscillator found in [I^. We show the 
contour lines of the Hamiltonian and the evolution of the 
uncertainty lines in Fig. |4^. 

The optimum orientation depends on the phase, i.e., 
Eq. (IMl) gives 

tan2()) = I f- - . (46) 

2 Va: Pj 

Expressing the phase space point {x,p) as x = A cos a, 
p = A sin a we get 

tan 2(/) = — cot 2a, (47) 

so that (j) = a ± 7r/4 which corresponds to the optimum 
orientation of the axes of the uncertainty ellipse inclined 
by ±7r/4 from the radius. 

The angular velocities of Eqs. (H51) and (H51) are 

= 8x{p‘^+x‘^), (48) 

Wc = 4x(/-ba:^). (49) 

Since Wc = uJv/2 and the state circles around the origin, 
(xrtPr) = (0,0), from Eq. we have 

{Xr,Pr) = l{Xo,Po). (50) 
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FIG. 5: Hamiltonian of the Jaynes Cummings type, Eq. (El) 
and motion of the uncertainty ellipse on the two branches. 


Thus, to keep the state close to the phase state point 
{xo,po) with the uncertainty ellipse optimally oriented, 
one needs to use the additional Hamiltonian 


Had = -4 :x{xI+pI) 


(51) 


/ xo\'^ ( PoV 

[x--) +(p-yj 


This Hamiltonian rotates the phase space around the 
point in the middle between the origin and the center 
of the uncertainty ellipse (a;o,Po) by twice the rate of the 
original rotation of (xo,po) around the origin. Contour 
lines of the resulting Hamiltonian H + H^d and evolution 
of the uncertainty lines are shown in Fig. |41 d. 


E. Jaynes-Cummings model 

Assume the Hamiltonian in the form 

= ( 52 ) 

The motivation comes from the Jaynes-Cummings model 
of a two-level atom interacting with a single mode field, 
the quantum Hamiltonian being 

Hjc = g{aa++a'a-) , (53) 

where g is a coupling constant and a± are the atomic 
raising and lowering operators. It was first shown in @ 
that this Hamiltonian can generate squeezed states of 
the optical field, which was elaborated in detail in [13- 
0. Let us assume the initial quantum state prepared 
as 

|$±) = |a) O(Ip) ±e*‘^|e)) (54) 

where |a) is the coherent state of the field with a ex¬ 
pressed as a = = 2~^l'^{x -P ip) and jp) and |e) 


(a) (b) 




FIG. 6: (Color online) Contour lines of the Jaynes-Cummings- 
like Hamiltonian and evolution of the uncertainty lines, (a) 
Hamiltonian ESI with the plus sign and with g = I, (b) the 
same Hamiltonian with the added term ESJ. The initial state 
is centered at (a;o,po) = (^i 0) and the uncertainty lines are 
plotted in times f = 0, 4, 8, and 16. 


are the ground and the excited atomic states, respec¬ 
tively. For times short compared to Trn/g the state re¬ 
mains approximately factorized so that one can study the 
evolution of the field separately from that of the atom. 
The mean energy of the state |<I>±) is (4>±|7?jc|$±) = 
ip-y/n = + ^ which corresponds to the 

classical Hamiltonian (EH)- 

Graph of the Hamiltonian (l5^ is a cone shown in 
Fig. m the two branches corresponding to the two signs 
of the atomic superposition in Eq. (IMl) . A state on 
the upper branch rotates clockwise whereas that on the 
lower branch counterclockwise. Note that in the quan¬ 
tum case, if the initial atomic state is different from 
2“^/^ (|p) ie*‘^|e)), the state evolves into a superposi¬ 
tion containing two separate coherent components of the 
field (i.e., a Schrodinger cat state, see [iM 1^). 

Each phase space point drifts along a circle of equal 
height. Since the magnitude of the cone slope does not 
depend on position (x,p), each point moves with equal 
speed = gj^/2. However, points on circles 

of different radii move with different angular velocities: 
points closer to the origin describe in the same time a 
bigger angle than points farther from the origin. As a 
result, a small area of the phase space is stretched in one 
direction and squeezed in the other one. We show the 
contour lines of Hamiltonian (15^ (branch with the plus 
sign) and evolution of the uncertainty lines in Fig. 

On using Eq. we find the squeezing rate as 


Q 


9 1 

\/2 yjx'^ +p'^ 


(55) 


Note that this rate agrees with the short-time value of 
squeezing evolution derived in the quantum model M- 
As can be seen, contrary to the Kerr model, the squeez¬ 
ing rate decreases with increasing the distance from the 
origin. Note also, that although the classical model is 
equally valid for any nonzero distance from the ori gin , 
the approximation derived in the quantum model 181 
























7 


works only for a;^ + ^ 1. 

The optimum orientation of the uncertainty ellipse 
given by Eq. (l24ll is 

tan2(;i = i--V (56) 

2 Va; pj 


A. Hamilton equations 

Assume Hamiltonian H{Jx, Jy, Jz) to be a differen¬ 
tiable function of Jx, Jy, and Jz- We postulate the 
Hamilton equations of motion as 


which is the same result as in the Kerr model. However, 
whereas in the Kerr model the distant part of the ellipse 
is ahead, in the Jaynes-Cummings model the part closer 
to the origin is ahead. 

For the angular velocities we find 


UJv = ± 


9 1 

•\/2 2 \ Jx '^ 


(57) 


and 


Wc = ± 


g 1 

\/2 +p'^' 


(58) 


with the upper (lower) sign corresponding to the upper 
(lower) branch of the Hamiltonian (1521) . Since uic = 2wt, 
and the state circles around the origin, {xh^Pr) = (0, 0), 
from Eq. dSH) we have 


{Xr,Pr) = {-Xo,-po), (59) 


J X 

OH 

dH 

(64) 

= 

dJy 



dH 

dH 

(65) 

Jy 

~ '^^Wz 

1 


dH 

dH 

(66) 

Jz 

~ 



These equations can be written in a condensed form as 


Ji — ^ijkJk 


dH 


(67) 


where eijk is the Levi-Civita symbol and Einstein sum¬ 
mation is used, or in the vector form as 


J = grad H X J. 


( 68 ) 


i.e., to keep the state close to the phase state point 
(xo,po) with the uncertainty ellipse optimally oriented, 
one needs to use the additional Hamiltonian 

Contour lines of the resulting Hamiltonian H + iJad and 
evolution of the uncertainty lines are shown in Fig. |B}d. 

VI. BLOCH SPHERE AS A PHASE SPACE AND 
SPIN SQUEEZING 

To describe dynamics of collective spin systems, one 
often depicts the states on a Bloch sphere with coordi¬ 
nates Jx, Jy, and Jz satisfying J^ + Jy + J^ = \J\^ where 
I J| is a constant. These numbers are related to angular 
momentum operators defined as 


^{d'^b + dP), 

(61) 

-^{d'^b- dP), 

(62) 

i(o'''a - Pb), 

(63) 


where a and b are the annihilation operators of two 
bosonic modes corresponding to the populations of atoms 
in two possible spin states. The angular momentum op¬ 
erators satisfy the relation J^ + Jy + J^ = y(y 
where N is the total number of particles. The Bloch 
sphere has properties of a compact phase space where 
the classical trajectories have been used, e.g., for Bohr- 
Sommerfeld quantization of spin Hamiltonians pH . [2H| . 
Here we use classical trajectories on the Bloch sphere to 
explore properties of various spin squeezing models. 


Equations (|64ll - (l66l) correspond in the quantum regime 
to the Heisenberg equations iA = [A, H] where combina¬ 
tions of operators Jx,y,z are taken in a symmetrical form 
and the commutation relations are {Jx,Jy\ = iJz with 
cyclical interchange of indexes. 

To see the correspondence of Eqs. (IMl) - (l66)l with the 


classical Hamilton equations in a planar phase space, let 
us assume N 1 and states near the pole with Jz~N/2 

and Jx,y N/2. 

Defining x = yj2/NJx, p = 

yYjNJy, 

and z = i/2/A J^ 

~ 12, we have 



dH lY dH 

(69) 

X 

dp N dz ’ 


dH lYdH 

(70) 

p 



For sufficiently weak dependence of H on J^ and large 
N, the second terms on the right hand side can be ne¬ 
glected and we arrive at the Hamilton equations for the 
ID motion. 


B. Integrals of motion and the Liouville theorem 

The equations of motion guarantee that both the 
Hamiltonian and -I- -I- are conserved quantities, 

dH/dt = 0 and d{Jx + Jy + Jl)/dt = 0. It is obvious es¬ 
pecially from Eq. (1551) : vector J moves perpendicularly 
to both the gradient of H and to itself. Thus, the phase 
space point moves on a surface of a sphere along lines of 
constant H with a speed proportional to the magnitude 
of grad H. 












One can also check that the Liouville theorem holds 
in this system. Let us assume a probability density 
p{Jx, Jy^ Jz,t) satisfying the continuity equation 


dp 


= -v-j, 


(71) 


where the current density is j = {Jx, Jy, Jz)p- Expressing 
the total time derivative of p as 


dp 

dt 


dp ■ dp ■ dp ■ dp 


using Eq. (ITT]) for dp/dt and Eqs. (I64l) - (l66l) for jx,y,z, 
one finds dp/dt = 0, i.e., the Liouville theorem holds. 
The probability density thus behaves as an incompress¬ 
ible liquid circulating along constant-Hamiltonian lines 
on a sphere. These results hint that the Bloch sphere 
with equations (l64ll ~ (l66ll represent a well-behaved phase 
space with classical evolution of state points. 


and 


^XX 

= V+ cos^ a -b 14 sin^ a. 

(76) 

^yy 

= V+ sin^ a -b 14 cos^ a. 

(77) 

Vxy 

14-1^- . „ 

= - - -Sin 2a, 

(78) 


where a is the orientation of the uncertainty ellipse. Ex¬ 
pressing the time derivatives by means of Eq. GU as 

Vxx = —‘^HzVxy + 2Jz{HxyVxx + HyyVxy), (79) 

Vyy = 2H2;Vxy - 2Jz{HxyVyy + HxxVxy), (80) 

- Hzi/^XX ^yy') ~b ^X^dlyy\/jy 77a;a;V^a;), (81) 

we find for the principal moments 

1^ = ^Jz [{ddyy ddxx) sin 2(x -j- 2Hxy cos 2o/ , 

(82) 


C. Evolution of moments 


which shows the dependence of the squeezing rate on the 
orientation. The optimum orientation occurs for 


Let us assume that the states are distributed in the 
vicinity of some phase-space point (ji°\ Jy^\ Ti°^) such 
that the values of Hamiltonian in any nearby point 
(-b AJx, + ^Jy, -b AJj) can be expressed 
by means of the Taylor expansion up to the quadratic 
terms. Denoting Jk = (Jk), Vm = {{Jk - Jk){Ji - Ji)), 
and the partial derivatives Hk = dH/dJk, etc., one finds 
(see Appendix [B] for details of the derivation) 


dJi 
dt 
dV^j 
dt 


^ijk (^dd^jJk “b Hjl^lk') ; 


Hi {^ilk^jk “b ^.jikVik^ 

“b diipJk {(lilk^jp “b ^jlk^ip) • 


(73) 


(74) 


Note that for the special case of quadratic Hamiltonians 
H = Wfe Jfc + XkiJkJi the results of Eo s (f7^ and (l7^ co¬ 
incide with the equations derived in for a quantum 
description of spin squeezing in Gaussian approximation, 
and in particular for Hamiltonian H = —ujJx -\- de¬ 
scribing a two-component Bose-Einstein condensate they 
coincide with the “Bogoliubov backreaction” equations of 

M- 


D. Squeezing rate and orientation of the 
uncertainty ellipse 


Let us first assume that the state is centered at the 
north pole of the Bloch sphere with Jx = Jy = 0, Jz > 0 
with no fluctuations in the radial direction, Vzk ~ 0, 
k = x,y, z. We express the variation matrix by means of 
the principal variances V±, where 


E± = 


14, 


^yy 


^{VxX-VyyY+Wjp, (75) 


tan2a= ^^^^^- , (83) 

for which we get t4i = ±Q14i with 

Q = \Jz\^{Hxx-HyyY+AHly, (84) 

which is analogous to Eq. (H. Note that if the coordi¬ 
nates are chosen such that Hxy = 0, the orientation of 
optimum squeezing corresponds to a = ±7r/4. 

In the case of a general position on the Bloch sphere 
one can proceed by first transforming the coordinate sys¬ 
tem to place the state to the pole and then use Eq. (1511) . 
Expressing the general result, one finds after some alge¬ 
bra 


Q = JTv{JH"JH" 


Tr 


(85) 


where H" is the matrix of the Hamiltonian second deriva¬ 
tives 


H" 


Hxx 

Hxy 

Hxz 

Hxy 

Hyy 

Hyz 

Hxz 

Hyz 

Hzz 


( 86 ) 


is the antisymmetric matrix corresponding to the co¬ 
ordinate vector J as 

/ 0 Jz —Jy \ 

J = \ -Jz 0 Jx (87) 

V Jy -Jx 0 j 


and I J| = Jx + Jy + Jz (s®® Appendix IQ for more de¬ 
tails of the derivation). 


2 
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E. Rotation of the Bloch sphere to keep the 
optimum squeezing orientation 


In general, the state driven by H is not only squeezed, 
but it also drifts along the Bloch sphere and gets the 
orientation of the uncertainty area rotated. Suppose that 
we want to keep the state at the chosen position and have 
its orientation optimal for fastest squeezing. Assuming 
that the initial distribution is sufficiently narrow that the 
terms HjiVik in Eq. d73l) can be neglected, compensation 
of the drift can be achieved by adding the Hamiltonian 


Hadi = —HkJk ( 88 ) 

= Wc • J, (89) 

where 

Wc = —gradi/, (90) 

and the derivatives of H are taken in {Jx,Jy,Jz) = 
{Jx: Jy^ Jz) • 

To find the rotation Hamiltonian that would keep the 
state optimally oriented, let us first consider that the 
chosen location is the north pole, i.e., Jx = Jy = 0, Jz > 
0 and the coordinate system is chosen such that H^y = 
0. This means that the optimum angle is a = 7r/4 which 
can be kept if variances of Jx and Jy are stretched with 
the same rate, i.e., we have V^x = Vyy in Eqs. GHD and 
o. This can be achieved if the Hamiltonian H + 
is supplemented with another term 


Had2 — i^vzJz 


(91) 


with 

^VZ = {Hxx + Hyy) . (92) 

This result is analogous to Eq. (1251) in the planar phase 
space. 

For a general position {Jx, Jy,Jz) on the Bloch sphere 
one can transform the coordinate system to get 

Had2 = ^ (^Hkk - jJi (93) 

= ujv ■ J, (94) 

where 

= (95) 

and the derivatives of H are taken in J = J. 



FIG. 7: (Color online) Phase space trajectories and evolution 
of the uncertainty region in the one-axis twisting scenario with 
Hamiltonian (iron . Quantities Jx,y,z are dimensionless here 
and in the next figures. 


corresponding to the one-axis twisting introduced in . 
Applying Eq. (IMI) we find 


Q 


2X 


Jx + Jy 

\j\ 


2x\J\siiJ0, 


(97) 


where Jx = |J|sin0cos0, Jy = |J|sin0sin0, and Jz = 
I J| cos 0. Thus, the fastest generation of squeezing occurs 
on the equator of the Bloch sphere with Q = 2x| J| corre¬ 
sponding to the quantum result Q = Nx with N = 2| J| 
being the total particle number. To compensate for the 
drift one rotates the Bloch sphere with Wc, 


We = 2x(0,0,J,), (98) 


and to keep the optimum orientation with 


(jj 


V 



Xsin^ 6J. 


J 


(99) 

( 100 ) 


Assume now an optimally located state at the equator, 
say, with Jy = Jz = 0 and J^, > 0. In this case there is no 
drift to compensate (wc = 0), and the optimal orientation 
of the uncertainty ellipse is kept by rotation with a;„ = 
x(|T|, 0, 0) so that the total Hamiltonian is 


F. Examples 


H,^, = H + Had 2 = xiJ! + \j\J.)- ( 101 ) 


1. One-axis twisting 

The simplest Hamiltonian used to generate spin 
squeezing is 

H = xJl (96) 


This Hamiltonian occurs in processes of classical bifur¬ 
cation studied, e.g., in [2^. We show the corresponding 
classical trajectories and the evolution of the uncertainty 
area in Fig. [71 These trajectories are equivalent to those 
found in the 6 — (j) plane by semiclassical analysis of the 
model in (2^. 
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FIG. 8: (Color online) Phase space trajectories and evolu¬ 
tion of the uncertainty region in the two-axis contertwisting 
scenario with Hamiltonian dMJ. 


FIG. 9: (Color online) Phase space trajectories and evolution 
of the uncertainty region in the Jaynes-Cummings model with 
Hamiltonian H + (wc-|-) ■ J of Eqs. (11081) , (IllOl) , and (nm, 
with the upper choice of sign. 


2. Two-axis countertwisting 


3. Spin squeezing by Jaynes-Cummings interaction 


Assume the Hamiltonian 


Let us assume the Hamiltonian 


H = x {Jl - 4) 


( 102 ) 


whose quantum counterpart corresponds to the two-axis 
countertwisting of 2^. Possible physical realization of 
such a Hamiltonian has been studied recently in 
For the squeezing rate we find 


Q = + (103) 

= xl J| \J sin^ 9 cos^ 2(j)-\- A cos^ 0, (104) 

which is maximized at the poles Jx = Jy = 0 with Q = 
2x1 J| and is zero at four points at the equator, = 0 
and (j) = ±7r/4, tt ± 7r/4. 

To compensate for the drift one needs to rotate the 
Bloch sphere with the angular velocity 


Wc = 2x(-J.,4,0), (105) 


and to keep the optimal orientation the sphere is rotated 
with 


Jx-Jy^ , ^ 

= -X |j|2 

= —X sin^ 0 cos 2())o^ (107) 

As can be seen, no rotation is necessary if the state is 
located at the optimum squeezing points Jx = Jy = 0, 
where Wc = = 0. This case is illustrated in Fig. [8] 

The results of the one-axis twisting and two-axis coun¬ 
tertwisting scenarios correspond exactly to the quantum 
results obtained from the Gaussian approximation in 
using the “twisting tensor” approach. 


H = ±g^J\J\-J,. 


(108) 


This form stems from the same considerations as in Sec. 
IV El assuming a two-mode field and an atom coupled 
to mode b of the field, the atom being prepared in pos¬ 
itive or negative superposition of the two levels. An¬ 
other physical realization proposed in is a collection 
of atoms with states |a) and |6), the latter being cou¬ 
pled by a laser field to a Rydberg state |r) for which the 
Rydberg blockade prohibits more than one atom in state 
|r). On using the definitions of x and p of Sec. IVI Al 
we see that near the north pole of the Bloch sphere with 
\Jx,y\ "C I J| = N/2, Jz ~ |J| Eq. (11081) reduces to Eq. 

dsil. 

For a general position, on using Eq. (1851) one finds 


g(| J| -I- Jz) 

4|J|V|J|-J,’ 


(109) 


which near the north pole of the Bloch sphere reduces 
to Eq. (1551) . For the drift compensation and for keeping 
optimum orientation of the uncertainty ellipse we find 
the angular velocities 


Wc 


UJy 


0,0,± r~-= 

\ 2./\J\^z, 

g I J| -I- Jz f 

■ 


( 110 ) 

( 111 ) 


We illustrate the classical trajectories and the evolution 
of the uncertainty area in Fig. [HI As can be checked, 
for states near the north pole of the Bloch sphere the 
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combined rotation is 

Wc + Wt, ~ — I _ ^ {~Jx, ~Jy} Jz) , (112) 

^8|J|(J2+J2) 

which means that the sphere is rotated around an axis 
intersecting the sphere oppositely to the state across the 
pole. This result is analogous to the rotation by Hamil¬ 
tonian dsni) in the planar phase space. 


VII. CONCLUSION 

The main result of this paper are equations m and 
(1851) giving the maximum squeezing rate Q in a planar 
phase space and on the Bloch sphere, respectively. In the 
planar phase space, Q is a function of second derivatives 
of the classical Hamiltonian which, in the zero-gradient 
points, is proportional to the difference of principal cur¬ 
vatures. On the Bloch sphere, formula (1551) gene ralizes 
the result found for quadratic Hamiltonians in where 
the maximum squeezing rate is proportional to the dif¬ 
ference of the maximum and minimum eigenvalues of the 
twisting tensor. The formulas with the second deriva¬ 
tives can be interpreted as using a local expansion of the 
Hamiltonian up to quadratic terms to “twist” the phase 
space neighborhood of the considered state. 

The other main results are the rotation frequencies of 
the phase space that keep the state at the right place 
and optimally oriented. They can be used as parame¬ 
ters of additional Hamiltonians to supplement the orig¬ 
inal Hamiltonian. These additional Hamiltonians them¬ 
selves do not produce squeezing and their addition does 
not influence the value of Q. They can be treated rather 
as instruments that optimize the exploitation, but do not 
change the amount of the resource. Their application 
transforms the point of interest into a saddle point with 
principal curvatures of equal magnitudes and opposite 
signs. 

It is interesting to note how several “quantum” re¬ 
sults could be found purely by classical means. Apart 
from the squeezing rates, these are, e.g., the “Bogoli- 
ubov backreaction” equations dZSl) and dZH) relevant for 
the description of two-component Bose-Einstein conden¬ 
sates. Although squeezing itself is sometimes described 
as a purely quantum phenomenon, we can see that it is 
not. What is “quantum” on squeezed states is rather the 
requirement on the minimum size of the uncertainty area, 
and sometimes the origin of the Hamiltonian governing 
the evolution (as, e.g., in the Jaynes-Cummings model). 
The classical results can be used for a quick estimation of 
the main properties of the states at the beginning of the 
squeezing process. At later stages of the evolution, how¬ 
ever, the quantum nature of our world starts revealing in 
the interference phenomena that cannot be described by 
the classical means. 



FIG. 10: Motion of the phase space point from {xi,pi) 
throngh {x 2 ,P 2 ) to ( 13 ,^ 3 ) expressed as rotation around 
{xr,Pr) by angle 7 . 
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Appendix A: Rotation radius and angular velocity of 
a phase space point 

Assume a phase space point moving from {xi,pi) 
through {x 2 ,P 2 ) to {xsjPs) as in Fig. [TUI First we ex¬ 
press {x 2 ,P 2 ) and (xsjPs) up to the second order of a 
short time interval dt. We have 


X2 ~ Xi + xidt +-xidt^ 



= xi + Hpdt+ ^(Ha;pHp - HppHa;)dt'^, 

(Al) 

P2 

1 2 

Ri Pi +Pidt + -Pldt 



= Pi- H^dt + - H^^Hp)dt^, 

(A2) 

where we 

have used 



~ ~ ^px^l HppPl 



— HxpHp HppHxf 

(A3) 


.. d 

Pi — — dixxdll HxpPl 

dt 



ddxxdip dixpdixt 

(A4) 

and the derivatives are taken in (xi,pi). Similarly we 

have 



X3 « 

1 2 

X2 + X2dt + -X2dt 

(A5) 

= 

X2 + HfUt + - Hj,fH(^'>)dt 

2 

1 

P3 ~ 

P2 + P2dt + ]^P2dt^ 

(A6) 


= P2- 
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(2) 

where the upper index in Hp , etc., means that the 
derivatives are taken in {x 2 ,P 2 )- We express these deriva¬ 
tives in terms of the derivatives in {xi,pi) with a correc¬ 
tion up to the first order in dt as 



— Hxpdp 



{.^xx^p ^xp^x^^^i 

(A7) 

H^p^ 

— Hp Hpxdx Hppdp 



— Hp “h {^HpxHp HppHx^dt^ 

(A 8 ) 


while for the second derivatives it is enough to keep only 

( 2 ) 

the zeroth order of dt, i.e., H^p = H^p, etc. We thus 
have 

a ^3 = Xi + 2Hpdt + 2{H,cpHp — HppHx)dt‘^, (A9) 
P3 = Pi - (AlO) 


A small angle 7 = uJcdt <C 1 can be expressed as 


7 = 


W X s 
H^Hpp 


- 2H^H„H. 


P 


p^-^xp 


m 




dt, (A13) 


from which we get Wc as in Eq. (E51) . 

Finally, we can express position of the center of the 
rotation as {xrjPu) = {xi,pi) + {Rx,Rp), where 


p _ 

^X - 5 

7 

(A14) 

n _ Sa: 

ttp — , 

7 

(A15) 


We can then express the vectors Si = ( 0 : 2 ,P 2 ) - {xi,pi) from which we get Eqs. ([26]) and (l27t . 
and S 2 = {x 3 ,p 3 ) — {x 2 ,P 2 ) of Fig. [TUlwhich are used to 
get their difference in = S 2 — as 


w — (^HxpRp RppRxiRxpR-x HxxRp)dt . (All) 

For the remaining calculations vectors 5 * 1^2 are taken up 
to the first order of dt as si ~ S 2 = s, i.e., 

s={Hp,-Hx)dt. (A12) 


Appendix B: Time derivatives of moments on the 
Bloch sphere 


Let us express the time derivative of Ji + AJi by means 
ofEq. dSll)' as 


J 


— {Ji -l- AJi) 
dt 


^ijk {Jk + AJfe) + HiAJi + —HisAJiAJg -I-... ^ 


— SijkJk + HjiAJi + —HjisAJiAJg -|-... 

d-tijk ^HjAJk + HjiAJ/^AJi + —HjisAJ/^AJiAJs -|-... ^ 


(Bl) 


If this expansion is taken in the mean value of {Jx, Jy, Jz) = {Jx, Jy, Jz), one can use (AJk) = 0 and (A A J/) = 14; 
so that the mean value of Eq. m yields 


— ^ijkdk (Rj ^^Rjls^ls ■ • 


+ Cfrfe ( HjiVkl + —HjisVkls + • . . 


(B2) 


where 14;s = (AJ^AJ^AJs), etc. If the distribution of {Jx, Jy, Jz) is sufficiently narrow around its mean, we can keep 
in each bracket just the first term so that we arrive at Eq. GSl). 

In a similar way we can express the time derivative of the product {Ji AJi){Jj -|- AJj) as 

^ [{Ji + AJi) {Jj -f AJj)] = [euk {Jj + AJj) + ejik {Ji + AJi)] {Jk + AJk) 

X -|- HipAJp —HiprAJpAJr + —id/prsAJpAA Jg -p ... ^ . (B3) 

Expressing the mean value of Eq. (IB3I) and using (IB2|) , one finds 


dV^j 

dt 


Rl {^ilk^jk “t” ^jlk^ik) “f RlpJk {^ilk^jp ^jlk^ip) 

~Hlpr {^ilk^jkpr ^jlk^ikpr) A ~^RlprsJk {^ilk^jprs A ^jlk^iprs) A . . . 


(B4) 
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For a sufficiently narrow distribution only the first two terms on the right hand side of Eq. (EH) are essential so we 
arrive at Eq. 


Appendix C: Squeezing rate for a general position on the Bloch sphere 

For a state centered at {Jx, Jy, Jz) the transformation 

( J'x\ 

= R2R1 



Jx 


Jy 

Jz 


(Cl) 


with 


R^ = 


y/ y/ 

~Jy _ 'R _ 

y/ y/ 

0 0 


0^ 

0 

1 / 


/ 


R2 = 


0 

\/ Jx^-\-Jy^ 


-y/jJ+jy \ 

yjy+jy^+R^ 

0 

_ R _ 

VTT+tiyjy J 


(C2) 


moves the state to the pole of the Bloch sphere with = J' = 0 where one can apply Eq. (IMl) as = 
|Jp [{Hx'x' — Hy'y') + AHliy,]. The derivatives with respect to the new coordinates Hx'x' = d^H/dJ'x etc are 
expressed using the chain rule 


etc, with 



( dJx 



fdJzY 


tc 

II 

tc 

) 

w) ’ 

dJx 

JxJz 

dJy 

JyJz 

dJx V Jx^ + Jy^ 

t>J’x 

y/Jx‘^+Jy'^y/Jx^+Jy^+Jz'‘ ’ 

iJJ'x 

y/Jx'‘+Jy'‘y/Jx'^+Jy‘^+Jz'^ ' 

9JL 

dJx 

— Jy 

C)Jy 

Jx 

— 0 
aj' 

^Jy 

y/Jx'^ + Jy'^ ’ 

dJ'y 

y/Jx'^ + Jy'^ ’ 

dJx 

Jx 

C)Jy 

Jy 

dJx _ Jx 


y/Jx'^ + Jy'^ + Jx'^ ’ 

dJ'^ 

^/Jx'‘ + Jy^ + Jx^ ’ 

^/Jx^+jV+H" 


Jfc we arrive at 


2 


On changing from the mean values to the coordinates of the point of interest Jk 
= \J\~"{{Jy+JzfHL + {Jx+J!yK + iJ^^ + JvTHl 

+4 (Jj + J2) (J2 + J2) + 4 (Jj + J2) (+ J2) Hi + 4 (+ J2) (+ Jl H, 

+2 \ JxJy ~ Jz {Jx + 'A/ + HxxHyy + 2 \j^ Jz ~ Jy {jX ^ ^ "^z)] HxxHzz 

+ 2 \^JyJz ~ Jx {Jx Jy + "^z)] JJyyJJzz 

-4.JxJy {Jy + j!) HxxHxy - UxJz {Jy + J^ HxxHxz + 4Jy J, (2 HxxHyz 

-AJxJy {Jl + Jl) HyyHxy + AJxJz ( + 2J^ + HyyHxZ " 4J, J, {j^ + J^ HyyHyz 

+AJxJy {Jl + Jl + 2Jl) HzzHxy - 4J,J, (jJ + J^) HzzHxz " AJyJz {Jl + Jl) HxzHyz 

-SJyJz {Jl + Jl) HxyHxz - SJxJz {Jl + Jl) HxyHyz - SJxJy {JI + JI) HxzHyz} , 

which can be abbreviated by using the Einstein summation as 

1 


= 


{{JsJs){^ijk^lpqJkJqHjlHpi) {^JkHy^qJq){^JsHgpjp) 


which corresponds to Eq. 


JwJw 

— 2{JkJk){JsHslHlqJq) + {JkJk){JsJs){HlqHql)\ 

— ^ijkJkHjlClpqJqHpi T — - (^(^ijkJk^jlpJpHlqCqysJs^rtvJpHti) • 
taking into account that Hij = ejjkJk- 


(C3) 


(C4) 


(C5) 


(C6) 


r 


[1] S. L. Braunstein, Phys. Rev. A 71, 055801 (2005). 


[2] J. Pefina, Quantum statistics of linear and nonlinear op- 















































14 


tical phenomena. Dordrecht Kluwer (1991). 

[3] R. Loudon and P. L. Knight, J. Mod. Opt. 34, 709-759 
(1987). 

[4] D. J. Wineland, J. J. Bollinger, W. M. Itano, and D. J. 
Heinzen, Phys. Rev. A 50, 67 (1994). 

[5] V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev. 
Lett. 96, 010401 (2006); Nature Photonics 5, 222 (2011). 

[6] P. Meystre and M. S. Zubairy, Phys. Lett. 89A, 390 
(1982). 

[7] G. J. Milburn, Phys. Rev. A 33, 674 (1986). 

[8] G. Drobny, A. Bandilla, and 1. Jex, Phys. Rev. A 55, 78 
(1997). 

[9] J. Bajer, O. Haderka, and J. Pefina, J. Opt. B: Quantum 
Semiclass. Opt. 1, 529 (1999). 

[10] J. Bajer and A. Miranowicz, J. Opt. B: Quantum Semi¬ 
class. Opt. 2, LIO (2000). 

[11] S. Wallentowitz and W. Vogel, Phys. Rev. A 53, 4528 
(1996); K. Banaszek and K. Wodkiewicz, Phys. Rev. 
Lett. 76, 4344 (1996); M. G. A. Paris, Phys. Lett. A 
217, 78 (1996). 

[12] G. Milburn and D. F. Walls, Opt. Gommun. 39, 401 
(1981). 

[13] G. M. Caves, Phys. Rev. D 23, 1693 (1981). 

[14] R. Tanas and S. Kielich, Opt. Gommun. 45, 351 (1983). 

[15] M. Kitagawa and Y. Yamamoto, Phys. Rev. A 34, 3974 
(1986). 


[16] J. Bajer, A. Miranowicz, and R. Tanas, Czech. J. Phys. 
52, 1313 (2002). 

[17] J. R. Kuklihski and J. L. Madajczyk, Phys. Rev. A 37, 
3175 (1988). 

[18] J. Gea-Banacloche, Phys. Rev. A 44, 5913-5931 (1991). 

[19] C. W. Woods and J. Gea-Banacloche, J. Mod. Opt. 40, 
2361-2379 (1993). 

[20] V. Buzek, H. Moya-Cessa, P. L. Knight, and S. J. D. 
Phoenix, Phys. Rev. A 45, 8190 (1992). 

[21] R. Shankar, Phys. Rev. Lett. 45, 1088 (1980). 

[22] A. Garg and M. Stone, Phys. Rev. Lett. 92, 010401 
(2004). 

[23] T. Opatrny, Phys. Rev. A 91, 053826 (2015). 

[24] A. Vardi and J. R. Anglin, Phys. Rev. Lett. 86, 568 

( 2001 ). 

[25] M. Kitagawa and M. Ueda, Phys. Rev. A 47, 5138 (1993). 

[26] T. Zibold, E. Nicklas, C. Gross, and M. K. Oberthaler 
Phys. Rev. Lett. 105, 204101 (2010). 

[27] A. Micheli, D. Jaksch, J. 1. Cirac, and P. Zoller, Phys. 
Rev. A 67, 013607 (2003). 

[28] T. Opatrny, M. Kolaf and K. K. Das, Phys. Rev. A 91, 
053612 (2015). 

[29] T. Opatrny and K. Mplmer, Phys. Rev. A 86, 023845 

( 2012 ). 



